<?php

set_time_limit(0);
class DifferentialEvolution {
	
	const D 			= 2;			// 2-dimensional array
	const CR 			= 0.3;			// Crossover probability; between 0..1
	const NP			= 1000;			// Population size; in many applications 10 x D
	const F				= 0.5;			// Mutation constant; between 0.5..1
	const TRIALS_COUNT 	= 50;
	
	private $data;
	private $population;
	
	public function __construct() {
		$this->data = array();
		$this->population = array();
	}
	
	public function perform() {
		$this->read('Task1.txt');
		
// 		$matrix = array(6, 7.5, 7.5, 13.75);
		
// 		$result = array(22.5, 41.25);
// 		$matrix = $this->inverse($matrix);
		
// 		var_dump(array($matrix[0] * $result[0] + $matrix[1] * $result[1], $matrix[2] * $result[0] + $matrix[3] * $result[1]));
// 		var_dump($this->method1());
		$result = $this->method2();
		var_dump($result[0]);
		var_dump($result[1]);
	}
	
	
	public function method1() {
		$n = count($this->data);
		$sumX = 0;
		$sumY = 0;
		$squareSumX = 0;
		$multXandY = 0;
		foreach ($this->data as $dataPoint) {
			list($x, $y) = $dataPoint;
			
			$sumX += $x;
			$sumY += $y;
			$squareSumX += pow($x, 2);
			$multXandY += ($x * $y);
		}
		
		$result = array($sumY, $multXandY);
		$matrix = $this->inverse(array($n, $sumX, $sumX, $squareSumX));
		
		list($c, $a) = array($matrix[0] * $result[0] + $matrix[1] * $result[1], $matrix[2] * $result[0] + $matrix[3] * $result[1]);
		var_dump($c); var_dump($a);
		
		return $this->calculateMSE($a, $c);
	}
	
	public function method2($meanSquareError = true) {
		$n = count($this->data);
		$sumX = 0;
		$sumY = 0;
		$squareSumX = 0;
		$multXandY = 0;
		foreach ($this->data as $dataPoint) {
			list($x, $y) = $dataPoint;
				
			$sumX += $x;
			$sumY += $y;
			$squareSumX += pow($x, 2);
			$multXandY += ($x * $y);
		}
		
		$this->generatePopulation();

		$count = self::TRIALS_COUNT;
		$best = array();
		$bestScore = 999999999;
		
		while ($count > 0) {
			foreach ($this->population as $i => $element) {

				list($a, $b, $c) = $this->generateCandidate($i);
				$trial = array();
				
				$k = rand(0, self::D-1);
				
				for ($j= 0; $j < self::D; $j++) {
					if ($this->getRandUni(0, 1) < self::CR || $j == (self::D - 1)) {
						// At least one vector from randomisation (mutation)
						
						$trial[$k] = $this->population[$c][$k] + self::F * ($this->population[$b][$k] - $this->population[$a][$k]);
					} else {
						$trial[$k] = $this->population[$i][$k];
					}

					$k = ($k + 1) % self::D;
				}
				
				
				$parameterA = $trial[0];
				$parameterC = $trial[1];
				
				if ($meanSquareError) {
					$score = $this->calculateMSE($parameterA, $parameterC);
				} else {
					$score = $this->calculateMedianSquareError($parameterA, $parameterC);
				}
				
				
				if ($score < $bestScore) {
					$bestScore = $score;
					$best = array($parameterA, $parameterC);
					$this->population[$i][$k] = $trial[$k];
				}
				
				unset($trial);
			}
			$count--;
		}
		
		return array($best, $bestScore);
	}
	
	public function generatePopulation() {
		for ($i = 0; $i < self::NP; $i++) {
			
			for ($j = 0; $j < self::D; $j++) {
				$this->population[$i][$j] = $this->getRandUni(-100, 100);
			}
		}
		
		return $this->population;
	}
	
	public function generateCandidate($parent) {
		$count = 0;
		$a = $this->getRandomPopulationInstanceIndex();
		while ($a == $parent) {
			$a = $this->getRandomPopulationInstanceIndex();
		}
		
		$b = $this->getRandomPopulationInstanceIndex();
		while ($b == $a || $b == $parent) {
			$b = $this->getRandomPopulationInstanceIndex();
		}
		
		$c = $this->getRandomPopulationInstanceIndex();
		while ($c == $b || $c == $a || $c == $parent) {
			$c = $this->getRandomPopulationInstanceIndex();
			
			$count++;
		}
		
		return array($a, $b, $c);
	}
	
	public function getRandomPopulationInstanceIndex() {
		
		return rand(0, count($this->population)-1);
	}
	
	/* Mean square error */
	public function calculateMSE($a, $c) {
		$n = count($this->data);
		$sum = 0;
		foreach ($this->data as $i => $element) {
			$d = $element[1] - ($a * $element[0] + $c);		// deviation
			$sum += pow($d, 2);
		}

		return ($sum / $n);
	}
	
	public function calculateMedianSquareError($a, $c) {
		$n = count($this->data);
		
		$deviations = array();
		foreach ($this->data as $i => $element) {
			$d = $element[1] - ($a * $element[0] + $c);		// deviation
			$deviations[] = pow($d, 2);
		}
		
		sort($deviations);
		if ($n % 2 == 0) {
			$x = (int)($n / 2);
			$median = ($deviations[$x] + $deviations[$x + 1]) / 2;
		} else {
			$x = ceil($n / 2);
			$median = $deviations[$x];
		}
		
		return $median;
	}
	
	/**
	 * 
	 * Input: 2x2 matrix as an 1-dimensional array
	 * 
	 * @param 	array $matrix
	 * @return 	integer
	 */
	private function detemrinant($matrix) {
		
		return ($matrix[0] * $matrix[3] - $matrix[1] * $matrix[2]);
	}
	
	private function inverse($matrix) {
		
		$determinant = $this->detemrinant($matrix);
		
		if ($determinant == 0) {
			throw new \InvalidArgumentException('Matrix determinant 0, cannot inverst matrix');
		}
		
		$newMatrix = array($matrix[3], (-1)* $matrix[2], (-1)*$matrix[1], $matrix[0]);
		foreach ($newMatrix as $i => $element) {
			$newMatrix[$i] = (1 / $determinant) * $element;
		}
		
		return $newMatrix;
	}
	
	public function read($fileName) {
		$file = fopen($fileName, 'r');
		
		$header = fgets($file);
		
		while (!feof($file)) {
			$line = fgets($file);			
			$row = preg_split('/\s+/', $line);
			
			if (count($row) < 2) {
				break;
			}
			
			$this->data[] = array($row[0], $row[1]);
		}
		
		fclose($file);
	}
	
	protected function getRandUni($min, $max) {
	
		return $min + mt_rand() / mt_getrandmax() * ($max - $min);
	}
	
}
?>


<?php 
$de = new DifferentialEvolution();
$de->perform();
?>